Naval  Oceanographic 

andAtmospheric  Technica;Note'‘ 

Research  Laboratory  March  1991 


ID 


These  working  papers  mere  prepared  forth*  finely 
dissemination  of  information;  this  document  does 
not  represent  the  official  pos&on  of  NQARL. 


ABSTRACT 


Archived  DMSP  Operational  Line  Scanner  data  are  stored  in  a  satellite  projection 
consisting  of  contiguous  blocks  of  line  and  pixel  information.  The  transformation  from  satellite 
line  and  pixel  coordinates  to  the  conventional  geodetic  latitude  and  longitude  is  derived.  A 
program  suit?*  ’e  for  use  on  a  personal  computer  is  provided.  The  program  can  be  used  for  fast 
registration  of  an  image  in  the  line  pixel  coordinate  system  to  any  other  projection. 
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Algorithms  for  Converting  Geodetic  Earth  Location 
to  Satellite  Time  and  Swath  Pixel  Coordinates 
for  the  DMSP  Satellite  System 

JL  Introduction 


The  purpose  of  this  note  is  to  document  fee  use  of  DMSP  fine  cphemeris  data  to  convert 
from  geodetic  latitude  and  longitude  to  an  Operational  Line  Scanner  (OLS)  time  and  swath  pixel 
number  suitable  for  navigation  and  mapping  to  an  arbitrary  projection.  The  DMSP  ephemeris 
is  generally  contained  in  the  DMSP  header,  and  is  available  with  all  data  supplied  by  AFGWC. 
The  technique  described  in  this  note  is  applicable  to  any  satellite  in  near  circular  orbit  with  a 
scan  pattern  perpendicular  to  the  cubital  subtrack. 

The  program  is  based  cn  the  original  routine  supplied  by  AFGWC  using  orbital  ephemeris 
data  from  the  Satellite  Data  Handling  System  (SDHS)  header  associated  with  the  satellite  data. 
The  orbital  elements  are  used  with  a  second  order  solution  to  Kepler’s  equations  to  find  the 
satellite  subtrack.  For  a  given  earth  location,  the  program  calculates  the  closest  great  circle 
distance  to  the  satellite  subtrack. 

A  C  program  for  the  implementation  of  this  code  is  provided  in  the  appendix. 

2.  Transformation  from  Lititude/Longitude  to  Line  Pixel 


2.1  Initialization 

The  transformation  from  geodetic  location  to  line  and  pixel  values  involves  initializing 
the  cubital  parameters  and  transforming  the  geodetic  point  to  a  satellite  coordinate  system.  The 
SDHS  header  provides  the  argument  of  the  epoch  latitude,  La ,  and  the  argument  of  perigee,  Lp. 
The  orbital  eccentric  anomaly  is  assumed  to  be  the  difference  of  the  latitude  and  perigee.  The 
initial  value  of  the  mean  anomaly  M,,  or  the  angle  moved  by  the  satellite  in  a  circular  oibit,  is 
evaluated  using  the  orbit  eccentricity,  e,  in  a  second  order  Kepler  equation. 


V  =L  -L 

a  a  p 

Mm  =  2e  sm^  -  -e2sin(2Fa) 
4 


(1) 


The  required  earth  location  latitude.  and  longitude,  are  transformed  to  a  cartesian 
coordinate  system,  with  Z  axis  pointing  to  the  North  Pole,  and  the  X  axis  at  the  ascending  node 
longitude  minus  the  required  longitude  (see  Figure  1).  The  transformation  is 
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Xr  =  cos(4>,)cos04) 


(la) 


yr  =  cos(4>/)sm(A) 


(lb) 


Zr  =  smft,) 


(lc) 


The  cartesian  coordinates  of  the  required  point  are  transformed  to  orbit  coordinates  (M.N)  by 
rotating  around  the  X  axis  through  die  inclination  angle  i,  (Figure  2) 


(2a) 

Mr  =  Frcos(0  +  Zrsin(tj 

(2b) 

This  completes  the  irit-alization  of  the  calculation.  The  following  sections  solve  for  the 
distance  of  closest  approach  between  the  satellite  subtract  and  die  required  point. 


2.2  Iteration  of  the  Satellite  Sabpoints 

The  smallest  angular  distance  between  the  required  point  and  the  satellite  subpoint  track 
is  found  by  a  projecting  of  the  required  point  into  the  orbital  plane,  and  then  iteratively  rotating 
the  satellite  in  its  orbit.  A  new  coordinate  system  fU,V>  is  defined  as  the  rotation  of  the  N  axis 
in  the  orbital  plane  by  an  angle  p.  The  coordinate  system  is  shown  in  Figure  3.  The  first  guess 
of  U  is  die  argument  of  latitude,  viz., 


Ur  =  Acos(p)  +  Mr  sin(p) 


(3a) 


!'r=-!Vrsm(p)+Mrcos(p) 


(3b) 


The  angle  between  the  U  axis  and  die  projection  of  the  required  vector  on  the  U  axis,  dp  is  given 
by 
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Figure  1.  Cartesian  orbit  coordinate  system. 
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d\ i  =  arctan 


(4) 


'Hi' 

Jr) 

Note  that  the  arc  tangent  in  the  above  equation  provides  an  angle  in  one  of  die  four  quadrants. 
Before  starting  the  iteration  to  minize  dn,  we  must  ensure  that  the  epoch  latitude  is  not  at  die 
projection  of  the  required  point,  i.e.,  dp  is  not  small.  If  this  is  the  case,  the  radial  distance 
between  the  satellite  subpoint  and  the  subtrack  is 

wr~-  Yrsin^  +  zrcos (0  { 5 ) 


and  we  do  not  need  to  iterate  farther.  Otherwise  we  iterate  the  satellite  vector  to  the  projection 
of  the  required  point  in  the  orbital  plane. 

The  iteration  proceeds  by  incrementing  the  angle  p  by  dp.  This  is  used  as  a  new  wtimatf 
of  the  eccentric  anomaly 


p7  =  p  -dp 

(6) 

(7) 

The  approximate  solution  to  Kepler’s  equation  is  given  by 

M  =  V -2esai(V)  -  -c2sin(2I0 


(8) 


from  which  the  satellite  angular  distance  increment  dM,  and  time  increment  dt  are 

dM  (9a) 

di  =dMINm  .  (9b) 


5 


where  N„m  is  the  satellite  mean  motion  in  radians  per  minute.  Daring  the  time  dt  the  earth  has 
rotated  an  angle  given  by 

\g  =  -dtA 


where  A  is  the  earth  rotation  rate. 


The  rotated  Cartesian  coordinate  system  is 

Xn  =  X^cosCAp  +  Tran(Ap  (Ha) 

rn  =  -  Xrsm(Ag)  +  Trcos(Ap  (1  «>) 

Zn  =  Zr.  (He) 

This  vector  is  now  rotated  by  the  inclination  angle  to  the  orbital  plane  (N,M,W  coordinates) 

N  -X  (12a) 

rt  n 

Mn  =  Yncos(i)  +  Znsm(i)  (12b) 

Wn  =  -  r„sm(i)  +  Z„cos(z) .  (12c) 


This  is  rotated  to  the  latest  guess  of  the  projection  of  the  required  point  on  the  orbital  plane  by 

un  =  tf„cos(p)  +  Mnsin(p)  (13a) 

V„  =  -  Nnd n(p)  ♦  Mncos{\i) .  (13b) 

From  this  we  calculate  the  latest  estimate  of  the  change  in  angle  from  the  current  U  axi~  to  the 
projection  of  the  required  point 
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dp 


=  arc  tan 


(14) 


If  this  angle  is  less  than  .0001  radians,  or  there  have  been  more  than  five  iterations,  the  loop  is 
terminated. 


23  Determination  of  Line  and  Pixel 

The  angular  distance  between  the  required  point  and  the  satellite  subpoint  is  given  by 

A=arcszn(Hg.  (15) 

In  case  the  loop  was  skipped  the  above  equation  is  evaluated  using  Wrt  from  equation  6. 

The  line  and  pixel  corresponding  to  the  required  earth  location  are  found  as  follow's. 
Assuming  that  line  number  at  the  argument  of  latitude  is  SLC,  the  scan  line  number  is  given  by 

SL  =  SL0  +  dt*SR  (16) 

where  SR  is  the  scan  rate.  The  pixel  number,  PK.  across  the  scan  is  given  by 

PIX  =  A  «  6  <17> 


where  8  is  satellite  polar  angle  between  pixels.  The  resulting  pixel  number  must  be  checked  to 
ensure  the  p&el  is  within  scan  limits. 

3.  Algorithm  Improvement 

In  general  this  approach  is  not  sufficiently  rigorous  to  provide  an  accuracy  of  1  pixel 
accuracy  (600  m)  in  earth  location.  The  general  approach  at  AFGWC  is  to  provide  two  sets  of 
third  order  polynomial  coefficients  determined  earlier  from  comparison  of  imagery  to  landmarks. 
This  possible  could  be  eliminated  by  using  a  more  correct,  iterative  approach  to  solving  the 
Kepler  equation  (Equation  9)  as  well  as  a  rigorous  determination  of  the  distance  from  the 
subtrack.  The  latter  calculation  should  use  higher  order  polynomials,  such  as  those  available  in 
SPG-4  (Hoots,  1980). 


7 


4. 


References 

Hoots,  and  R.L.  Roerich,  1980:  Models  for  Propagation  of  NORAD  Element  Sets,  Project 
Space  Track  Rpt.  3,  Aerospace  Defense  Command,  Peterson  AFFB,  Colorado  Springs,  CO  80840. 


APPENDIX  A 

Turbo  C  Code  for  Converting  Geodetic  Coordinates  to  Time  Pixel  Coordinates 


/*  LP2LL 

11-5-89 
A.  Goroch 

Convert  Latitude,  Longitude  to  Line  Pixel  for  DMSP  Fine 
Use  AFGWC  standard  code 

*/ 


#include  <stdio  h> 

#include  <malh.h> 

#define  PI  3.14159265358 
#deflne  D2R(A)  (A*PI/180.0) 

#define  PIXELPFRRADIAN  8*22*258 
#deflne  SCANLINEPERMINUTE  704 

typedef  struct  { 

float  ArgLatitude; 
float  AigPerigee; 
float  Eccentricity; 
float  LAN; 
float  Inclination; 
float  N; 

float  RelEaithRR; 

}  DMSPFINEORBIT; 
typedef  struct  { 
int  Pixel; 
int  Line; 
float  Latitude; 
float  Longitude; 

}  SATGEOLOCATION; 

void  CnvtGeodetic2PolOib(  DMSPFINEORBIT  Orbit,  SATGEOLOCATION  *SGL); 
main() 

I 

I*  Driver  routine  to  test  out  location  calculation  */ 

/*  Initialize  parameters  */ 


f*  Argument  of  latitude  in  Degrees  */ 
l*  AOP  in  degrees  */ 

/*  no  units  */ 

f*  Longitude  of  Ascending  Node  (Degrees  )*/ 
l*  inclination  (degrees)*/ 
l*  Anomalistic  mean  motion  */ 

/*  Relative  Earth  rotational  rate  */ 


I*  Degrees  */ 

/*  Degrees  East  */ 
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DMSPHNEORBIT  Oib  =  {12.0, 270.0,  .000155, 270.0,  99-15, 14.127,  -25.5}; 
SATGEOLOCATION  Sgl; 
float  Lat,  Lon; 

for(  Lat  =  -80.0  ;  Lat  <  80.0  ;  Lat  +=  10.0  ) 

{ 

over  latitudes  */ 

for(  Lon  =  0.0  ;  Lon  <=  360.0  ;  Lon  +=  10.0  ) 

{ 

Longitudes  */ 

Sgl. Latitude  =  Lat; 

Sgl.Longitude  =  Lon; 

CnvtGeodetic2PolOrb(  Orb,  &Sgl); 
printff  Lat  %5.1f.  Lon  %5.1f,  Line  %7d  Pixel  %7d\n". 

Sgl -Latitude,  Sgl.Longitnde,  Sgl-Line.  Sgl -Pixel); 

[*  print  out  LatLon,  Line  Pixel  */ 

}  /*End  of  longitudes  */ 

}  [*  end  of  latitudes  */ 

}  [*  End  of  main  */ 

I* _ CnvtGeodetic2PolOrb _ *1 

void  CnvtGeodetic2PolOrb(  DMSPHNEORBIT  O,  SATGEOLOCATION  *SGL) 

{ 

float  V0,V, 

M0,  M,  DM, 

SL  Ci, 

Ecc_2,  Ecc_sqr_75, 

A, 

X.  Y,  Z, 

RM,  RN,  RW, 

U,  DU,  DT, 

RU,  RV, 

Xt,  Yt,  Zt, 

DistSubTrack  /* Angular  distance  to  subtrack  */ 

* 

int  I  =  0,  Numlterat:ui»3  =  5  'J*  Index  and  max  no  of  iterations  */ 


V0  =D2R(0  AigLatitude)  -  D2R(0  ArgPerigee); 

Ecc_2  =  2  *  O.Ecccntricit>-; 

Ecc_sqr_75  =  Ecc_2  *  3  /  8  *  O -Eccentricity; 

M0  =  V0  -  Ecc_2  *  sin(  V0  )  +  Ecc_sqr_75  *  sin(  2  *  V0  ); 
Si  =  sin(  D2R(  O -Inclination  )  ); 

Ci  =  cos(  D2R(  O .Inclination  ) ); 


f*  loop 
l*  Loop  over 
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[*  preliminaries  completed  find  the  closest  line,  pixel  */ 


f*  Locate  Requested  vector  in  orbital  plane  wit  ascending  node  */ 

A  =  D2R(  O.LAN  )  -  D2R(  SGL->Longitude  ); 

X  =  cos(  D2R(  SGL->Latitude  ))  *  cos(A); 

Y  =  cos(  D2R(  SGL->Latitude  ))  *  sin(A); 

Z  =  sin(  D2R(  SGL->Lathude  )); 

f*  rotate  angle  i  around  X  to  get  to  orbital  coordinates  */ 

RN  =  X: 

RM  =  Y  *  Ci  +  Z  *  Si; 

l*  Get  initial  estimate  for  angle  satellite  vector  sut  move  to  coincide  with 
projection  of  R  on  UW  plane  */ 

U  =  D2R(O.ArgLalhude)  ; 

RU  =  RN  *  cos(  U  )  +  RM  *  sin(  U  ); 

RV  =  -RN  *  sin(  U  )  +  RM  *  cos(  U  ); 

DU  =  atan2(RV,  RU): 

DT  =  0.; 

1  =  0; 

RW  =  -Y  *  Si  +  Z  *  Ci; 

while{  (  fabs(  DU  )  <  1.0e-4  )  &&  ( I  <  Numlterations  )) 

{  f*  loop  until  convergence  or  four  tries  *1 
U  +=  DU: 

V  =  U  -  OArgPerigee; 

M  =  V  -  Ecc_2  *  sin(  V  )  +  Ecc_sqr_75  *  sin(  V  *  2  ); 

DM  =  M  -  MO; 

DT  =  DM  /  O  N; 

J*  rotate  (X,  Y,  Z)  around  Z  through  change  in  longitude  (  Dt  *  RERR)  */ 


A  =  -  DT  *  O-RelEarthRR; 

Xt  =  X  *  co s(  A  )  +  Y  *  sin(  A  ): 
Yt  =  -  X  *  sin(  A  )  +  Y  *  cos(  A  i: 
Zt  =  Z; 
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f*  Rotate  new  location  around  X  through  I  */ 


RN=  Xt; 

RM  =  Yt  *  Ci  +  Zt  *  Si; 

RW  =  -Yt  *  Si  +  Zt  *  Ci; 

I*  Rotate  through  satellite  rotation  */ 

RU=  RN*cos(U)  +  RM*sin(U); 

RV  =  -RN  *  sin(  U  )  +  RM  *  cos(  U  ); 

I*  Get  the  new  DU.  angle  from  U  to  R  on  UV  plane  */ 

DU  =  atan2(RV,  RU); 

1++;  /*  Increment  the  counter  */ 

}  I*  end  of  iteration  over  oibh  locations  *f 

DistSubTrack  =  asin(  RW  ):  l*  Angular  distance  from  subtrack  */ 

printff  Si  %f  Ci  %f,  NnR  =(  %f,  %f,  %f)\n  Rt=(  %f.  %f.  %f)\t  RW=  %f  subtrack  %f 
Sn", 

Si,  Ci, 

X,  Y,  Z.  Xt.  Yt,  Zt,  RW,  DistSubTrack); 

f*  Convert  Radial  distance  and  time  increment  to  pixel  scan  line  */ 

SGL->Pixel  =  DistSubTrack  *  PIXELPERRADIAN  ; 

SGL->Line  =  DT  *  S  CANLINEPERMINUTE  ; 

}  /*  end  of  CnvtGeodetic2PoIOib  */ 
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APPENDIX  B 


DMSP  Constants  of  Motion 


Hie  following  constants  are  assumed  to  characterize  the  physical  characteristics  of  the 
OLS  sensor  in  fine  mode. 

The  satellite  polar  angle  (5)  between  adjacent  pixels  is  converted  to  the  central  angle 
(earth  based)  subtended  by  the  sensor 


ptrH 


f  =arc  sin 


where  h^,  is  the  satellite  altitude  above  the  mean  earth  surface  in  m,  and  is  the  mean  earth 
radius  in  m.  The  DMSP  scan  rate  is  given  as  704  scans  per  minute. 
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